# -*- coding: utf-8 -*-
"""
#-------------------------------------------------------------------------------
# Nom:         optimGR4J
#
# Auteur:      Alain Gauthier
#
# Créé le:     23/08/2013
# Copyright:   (c) eaucéa 2014
# Licence:     GNU Public Licence (GPL)

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

Module permettant l'optimisation des paramètres de GR4J selon des critères
spécifiés ('Nash' ou bien critère régional)

#-------------------------------------------------------------------------------
"""

import datetime
import numpy as np
import numpy.linalg as LA
import scipy.optimize as opt
import analyseChroniques
import typeCritereNash
import gr4j

#-------------------------------------------------------------------------------

def lancerOptimisation(AcritereEfficacite,
                       AchroniquesPluie,
                       AchroniquesEvapo,
                       AchroniquesQmesure,
                       AdonneesStation,
                       fichierLog):
    """
    Optimisation des paramètres de GR4J pour la station de mesure.
    AcritereEfficacite : identifiant de la fonction Nash à optimiser
    AchroniquesPluie : chronique de pluie sur le BV de la station
    AchroniquesEvapo : chronique d'evapotranspiration sur le BV de la station
    AchroniquesQmesure : chronique des débits mesurés par la station de mesure
    AdonneesStation : dict de paramètres nécessaires au calage (Xparam : valeurs initiales des paramètres X de GR4J,
    duree_eval : (debut,fin) intervalle de temps de l'évaluation du critère d'efficacité,
    aire : aire du BV au niveau de la station)
    """

    # initialisations

    # paramètres à optimiser et paramètres fixés
    xparam = [int(x) for x in AdonneesStation['Xparam']]
    Lxini=np.array([xparam[0], xparam[2], xparam[3]]) # X1, X3, X4 à optimiser
    LXautre = np.array([xparam[1]]) # X2 valeur fixée

    # niveau initial de remplissage des réservoirs S/X1 et R/X3
    LniveauIni = np.array([0.6,0.7])

    dt_station = [analyseChroniques.decoder_date(tps) for tps in AdonneesStation['duree_eval']]

    # chronique de débit mesuré (en mm)
    tab_dtPluie = AchroniquesPluie['temps']
    dateP0 = tab_dtPluie[0]
    dateP1 = tab_dtPluie[1]
    dtPluie = dateP1-dateP0
    # cas particulier du pas de temps mensuel : fixer dtPluie=31 jours pour la fonction de conversion
    if dateP1.month != dateP0.month:
        dtPluie = datetime.timedelta(31)

    LQmes = analyseChroniques.debitEn_mm(AchroniquesQmesure,
                                         dtPluie,
                                         AdonneesStation['aire'])

    # alignement du début des chroniques :
    LchroniqueEvapo,LchroniquePluie,LQmes = analyseChroniques.alignerChroniques(AchroniquesEvapo,AchroniquesPluie,LQmes)

    # chronique de pluie (en mm)
    LPluie = np.array(LchroniquePluie['valeur'])
    # chronique d'évapotranspiration (en mm)
    Levapo = np.array(LchroniqueEvapo['valeur'])

    # chronique des débits mesurés en ayant retiré les valeurs négatives pour la simulation
    LQ_calage = np.array([np.max([val,0]) for val in LQmes['valeur']])
    # intervalle de temps d'évaluation du critère d'efficacité (indice dans le tableau des chroniques de pluie)
    tab_dtPluie = LchroniquePluie['temps']
    Ldt_eval = np.array([tab_dtPluie.index(dt_station[0]),tab_dtPluie.index(dt_station[1])])

    # test de validité de l'intervalle d'évaluation
    if np.min(LQ_calage[Ldt_eval[0]:Ldt_eval[1]])<0:
        raise ValueError('Tentative de calage sur un intervalle manquant de donnees : '+str(dt_station))
    # optimisation
    LniveauIni = np.array([0.6,0.7])

    fichierLog.write('param à optimiser : ' + str(Lxini)+'\r\n')
    fichierLog.write('dt_station = '+str(dt_station)+'\r\n')
    fichierLog.write('Ldt_eval = ' + str(Ldt_eval)+'\r\n')
    fichierLog.write('tps_eval = ' + str(LQmes['temps'][Ldt_eval[0]])+ ' , ' + str(LQmes['temps'][Ldt_eval[1]-1])+'\r\n')
    fichierLog.write('mean(LQ_calage) = ' + str(np.mean(LQ_calage[Ldt_eval[0]:Ldt_eval[1]]))+'\r\n')
    fichierLog.write('min(LQ_calage) = ' + str(np.min(LQ_calage[Ldt_eval[0]:Ldt_eval[1]]))+'\r\n')
    fichierLog.write('max(LQ_calage) = ' + str(np.max(LQ_calage[Ldt_eval[0]:Ldt_eval[1]]))+'\r\n')
    LXfmin = opt.fmin(fmin_gr4j_station,
                      Lxini,
                      (LXautre,LniveauIni,Ldt_eval,LPluie,Levapo,LQ_calage,dtPluie,
                      AcritereEfficacite,fichierLog),
                      maxiter=1000)

    # résultat
    resultat = [LXfmin[0],LXautre[0],LXfmin[1],LXfmin[2]]

    fichierLog.write('optimisation terminée\r\n')
    fichierLog.write('X optimal = ' + str(resultat)+'\r\n')

    # paramètres GR4J optimisés et fixés
    return resultat

#-------------------------------------------------------------------------------

def fmin_gr4j_station(x,AXautre,AniveauIni,Adt_eval,APluie,Aevapo,AQ_calage,AdtPluie,AcritereEfficacite,fichierLog):
    """
    fonction de calcul de l'erreur du débit par GR4J par rapport aux données Qj
    lors d'une simulation
    Les paramètre de GR4J sont répartis de la manière suivante:
    x[0], x[1], x[2] # X1, X3, X4 à optimiser
    AXautre = valeur # X2 valeur fixée

    paramètres :
    x          = variable à optimiser : paramètres de GR4J sauf AXautre
    AXautre    = paramètres de GR4J de simulation
    AniveauIni = vecteur contenant les niveaux de remplissage initiaux des réservoirs (en %) de GR4J
    Adt_eval   = intervalle de temps de simulation pour le calcul de l'erreur
    APluie     = chronique de pluie
    Aevapo     = chronique d'évapotranspiration
    AQ_calage  = chronique de débit mesuré pour calcul d'erreur
    AdtPluie   = pas de temps des chroniques (en jours)
    AcritereEfficacite = critère d'éfficacité choisi pour le calcul de l'erreur
    """

    # initialisation pour simulation GR4J
    lX=np.array([x[0],AXautre[0],x[1],x[2]])
    # pas de temps en heures
    Lpas=24*AdtPluie.days

    fichierLog.write('dans fmin\r\n')

    fichierLog.write('lX = ' + str(lX)+'\r\n')
##    fichierLog.write('APluie = '+str(APluie[:10])+'\r\n')
##    fichierLog.write('Aevapo = '+str(Aevapo[:10])+'\r\n')
##    fichierLog.write('AniveauIni = '+str(AniveauIni)+'\r\n')
##    fichierLog.write('Lpas = '+str(Lpas)+'\r\n')

    # appel de GR4J
    resultatGR4J = gr4j.gr4j(lX,APluie,Aevapo,AniveauIni,Lpas)

    # calcul de l'erreur selon le critère d'efficacité AcritereEfficacite
    # pas de transformation
    LQ_GR4J = resultatGR4J['Q'][Adt_eval[0]:Adt_eval[1]]
    LQ_ref = AQ_calage[Adt_eval[0]:Adt_eval[1]]

##    fichierLog.write('LQ_ref = ' + str(LQ_ref[:10])+'\r\n')
##    fichierLog.write('mean(LQ_ref) = ' + str(np.mean(LQ_ref))+'\r\n')
##    fichierLog.write('min(LQ_ref) = ' + str(np.min(LQ_ref))+'\r\n')
##    fichierLog.write('max(LQ_ref) = ' + str(np.max(LQ_ref))+'\r\n')
    Nash = calculerNash(LQ_GR4J,LQ_ref,AcritereEfficacite)

    erreur = 1-Nash

    fichierLog.write('Nash = '+str(Nash)+'\r\n')
    fichierLog.write('erreur = '+str(erreur)+'\r\n')

    # vérification des bornes de x
    if (np.min(x)<0):
        erreur = 5

    # résultat
    return erreur

#-------------------------------------------------------------------------------

def calculerNash(AQ_GR,AQ_ref,AcritereEfficacite):
    """
    calcul du critère de Nash selon le paramètre AcritereEfficacite.
    AQ_GR est la chronique de débit simulé
    AQ_ref est celle de débit mesuré
    """
    # exclusion des valeurs de AQ_ref (les simulations ne génèrent pas de valeurs nulles)
    id_ok=np.where(AQ_ref>0)
    # Q
    if len(id_ok[0]>0):
        lQ_GR = np.array(AQ_GR[id_ok[0]])
        lQ_ref = np.array(AQ_ref[id_ok[0]])
    else:
        lQ_GR = np.array(AQ_GR)
        lQ_ref = np.array(AQ_ref)
    # sqrt(Q)
    if AcritereEfficacite == typeCritereNash.listeCritereNash.index(typeCritereNash.__NASHVQ__):
        lQ_GR = np.sqrt(lQ_GR)
        lQ_ref = np.sqrt(lQ_ref)
    # log(Q)
    elif AcritereEfficacite == typeCritereNash.listeCritereNash.index(typeCritereNash.__NASHLOGQ__):
        lQ_GR = np.log(lQ_GR)
        lQ_ref = np.log(lQ_ref)

    # calcul du critère de Nash
    QrefMoyen=np.mean(lQ_ref)
    Nash=1-(LA.norm(lQ_GR-lQ_ref,2)**2)/(LA.norm(lQ_ref-QrefMoyen,2)**2)
    # résultat
    return Nash

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
